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Abstract 

Tomographic techniques are successfully applied to obtain 4D images of the 
tropospheric refractivity in a local dense network. In the lower atmosphere both 
the small height and time scales and the non-dispersive nature of tropospheric 
delays require a more careful analysis of the data. We show how GPS data is 
processed to obtain the tropospheric slant delays using the GIPSY-OASIS II 
software and define the concept of pseudo-wet delays, which will be the observ- 
ables in the tomographic software. We then discuss the inverse problem in the 
3D stochastic tomography, using simulated refractivity fields to test the system 
and the impact of noise. Finally, we use data from the Kilauea network in Hawaii 
and a local 4x4x41-voxel grid on a region of 400 Km 2 and 15 Km in height to 
produce 4D refractivity fields. Results are compared with ECMWF forecast. 



1. Introduction 

One goal of ongoing research on the application of tomographic techniques to 
the atmosphere using GPS signals is to obtain 4D images of the refractivity in the 
troposphere using a ground network. In this paper we present our first results on 
the subject, showing that a careful combination of GPS data processing together 
with tomographic Kalman-filtering can reproduce the state of the refractivity in 
the neutral atmosphere. 

The effect of the atmosphere on GPS signals is measured as an extra delay. 
The ionospheric electron content produces a dispersive delay and can hence be 
corrected by using two different frequencies and linearly combining them. The 
neutral atmosphere, however, induces a delay independent of frequency. This 
delay (AL) can be related to the refractivity by 

AL = J lO^Ndl, (1.1) 

JVw 77.6^ + 3.73- 10 5 ^, (1.2) 
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where P is the total atmospheric pressure (mbar), P w is the water vapor pressure 
(mbar), and T is the atmospheric temperature (K) ( \Smith et al. 1953| ). The 
extraction of this slant delay from GPS measurements requires accurate mod- 
eling. To take into account the dependence of the tropospheric delay on the 
satellite elevation, mapping functions are used, and a Zenith Total Delay (ZTD) 
is estimated at each station. To consider azimuthal variability, a horizontal gra- 
dient is also estimated. The time evolutions of these parameters is modeled as 
random walk stochastic processes. The total atmospheric delay can be parti- 
tioned into two components: the hydrostatic delay, due to the dry gases in the 
troposphere and the nondipole component of water vapor, and the wet delay, 
due to the dipole component of water vapor. The contribution of the hydrostatic 
component to the ZTD is larger than that of the wet and it can be estimated from 
surface pressure measurements. We remove the zenith hydrostatic delay (ZHD) 
from the total slant delays through an accurate estimation of the pressure ac- 
cording to site measurements, and then form the pseudo-wet delays (PWD): we 
map the time series zenith wet delay (ZWD) and total horizontal gradients back 
to each ray direction and add the residuals to account for any data mismodeling. 
The refractivity associated is termed pseudo-wet refractivity (N). The PWD's 
are the observables in the tomographic approach. The reason for removing of 
the ZHD is two-fold. On the one hand, it reduces the scale height of the problem 
(typically, the hydrostatic component has a scale height of 10 Km while the wet 
component has a 2 Km scale height) and thus eases the numerical solution. In 
addition, the N field obtained is closely related to the water vapor distribution. 

In the paper we will first discuss the approach taken to estimate the PWD's 
using the JPL software GIPSY-OASIS II (GOA II) (see \\Webb et al. 199% ) 
and the network used in the study. Then we briefly discuss the inverse problem 
and analyze the impact of the needed constraints through simulations. Finally, 
tomographic solutions with GPS data measured in the Kilauea network in Hawaii 
during 1st February 1997 and their comparison with independent measurement 
sources will be presented. 



2. Tropospheric Slant Delay Estimation 

The GPS signals L\ and L2 are first combined to remove the effect of the 
ionosphere. The resulting observable, Lq, is modeled in GOA II using a Point 
Positioning Strategy in terms of tropospheric effects (ZWD+ZHD+Gradients), 
geometrical factors and receiver clocks. Satellite clock corrections and orbits are 
provided by JPL, as well as earth orientation parameters. Tropospheric delays 
are mapped to the zenith using Niell mapping functions ( ^Niell 1996 1) and related 



to the gradients following the expression ( \\Davis et aZ.T993| ): 

Dl = m h (e)D hz + m w (e)D wz + 

+ mA(e) cot e[GN cos <fi + Ge sin<^], (2-1) 

The time dependence of the tropospheric component is handled using a kalman 
filter and a random walk stochastic process, with a drift rate of 0.25 cm/\/7i for 
the zenith delay and 0.03 cm/\//i ( [ Bar-Sever et al. 1998|] ) for the gradients. 



We then form the PWD from the solution as discussed above. 



2.1. Network and data considered 



GPS data was tracked with a sampling of 30 s in the Kilauea network in Hawaii. 
A map with the stations is shown in Figure |3[ It can be seen that they cover an 
area of about 400 Km 2 . The heights of the stations are distributed in a range of 
about 2000 m, which seems particularly well-suited for vertical resolution in the 
tomographic solution. 



3. Tomographic approach and results 
3.1. The inverse problem 

In tomography, one wants to obtain the solution fields (refractivity in the tro- 
pospheric case) from the integrated value along a ray path. The amount and 
distribution of data, however, is often insufficient to determine a unique solu- 
tion. Thus, some additional information has to be added to the system. To this 
end, we rewrite the pseudo-wet refractivity as N(f,t) = J2j ajit)^ jif) + e if^) 
where \J/j(r) is any set of basis functions and e(f,t) is the quantization error. In 
our tomographic approach, three dimensional voxels are used. Then, our obser- 
vations are modeled by %n = Tj = J sl N(f, t)dl = J2j a j(t) I s i ^j{f)dl + £(r, t). 
This is a set of linear equations of the form Ax = y. As mentioned, the sys- 
tem may not have a solution and thus, we seek to minimize the functional 
X 2 (x) = (y — Ax) T ■ (y — Ax). Even so, the unknowns may not be fixed by 
the data, although the number of equations is usually much greater than the 
number of unknowns. This is so because we are trying to obtain a solution with 
more degrees of freedom (more voxels) than the resolution carried by the data. 
In previous work, (see [Ruffini et al. 1998 1 and \Rius et al. 1997 |) we have 



discussed the use of a correlation functional to confine the spatial spectrum of 
the solution to the low portion of the frequency space. The same concept can 
be expressed by adding new equations (constraints) that impose that the density 
in a voxel be a weighted average of its neighbours. Let us now take a closer 
look to the constraints. Horizontal smoothing is needed to account for the non- 
uniform distribution of the rays. Water vapor is mostly concentrated in the first 
2 Km in height. Now, if we aim for a tomographic solution with a resolution of 
350 m (a trade-off between data resolution and computer load), some smooth- 
ing constraints have to be added in the vertical direction. Finally, we impose 
a zero value constraint in the highest layer. The constraints are represented by 
the equation B ■ x = 0. To smooth out the time variability, a Kalman filtering 
was used, modeling the troposphere as a random walk stochastic process with a 
drift rate of 5 = 0.14 cm/(Km yh). The tomographic process is integrated in 
our home-made software package LOTTOS (LOcal Tropospheric TOmography 
Software). 



3.2. Results with simulated data 

In order to test the software for the network's geometry and tune the different 
parameters, simulations play an essential role. We have simulated a 3D pseudo- 



wet refractivity field following the expression (following \Davis et al. 1993| ]) 

N(p, z) = N^e~^(l +g w -p)+ iVoVMl + 9d • p) (3.1) 

where h w = 2 Km, hd = 10 Km, and g w are the horizontal gradients up to a 
multiplicative constant in the hydrostatic and wet component, respectively and 
p is the horizontal coordinate. We have set = 150 mm/Km and N^. = 2 
mm/Km, the latter to account for any residual hydrostatic component, gj = 
and two different non-zero g w values, applied depending on the latitude of the 
station. The geometry of the rays for the 1st February 1997 in the Kilauea 
network has been used to generate the simulated rays. The PWD for each ray 
has been formed according to L s { m = J2vox Ni ■ k, where iVj is the pseudo-wet 
refractivity associated to voxel i and li is the length of the ray across this voxel. 
The weight of the constraints was tuned to give the best fit with simulations, 
successively adding the constraints in the horizontal plane, then the vertical 
smoothing constraints, and finally adding a zero value constraint for the highest 
layer. In Figure |] we show the error in reconstructing the simulated field. 

Vertical reconstruction has also been verified adding a perturbation at dif- 



ferent heights to the simulated fields of Equation 3.1; good vertical resolution 
is achieved thanks to the distribution in height of the stations. We also note 
(see Figure [l]) that vertical reconstruction is gradually lost as the perturbation 
is placed at higher altitudes (above 2 Km, resolution starts degrading). How- 
ever, when the boundary constraint kicks in, the solution agrees again with the 
simulated field, slowly decaying to zero. 

The impact of noise in the rays has also been considered. Following the data 
processing described above, we have seen that post-fit residuals are dependent 
on elevation (roughly as l/sin(e), where e is the elevation); therefore, we can 
assume a noise value for the ZWD, then map it to the slant directions, and add 
it to the total simulated delay. We show in Figure ||| the results of the noise 
analysis. As it can be seen, a ZWD noise of 0.5 cm represents a noise in the 
solution of less than 3 mm/Km rms (about 2% of the maximum value). 



3.3. Results from the Kilauea network 

Pseudo-wet delays from the Kilauea network in Hawaii have been input to the 
system. We have processed data from 18 stations for February 1st, 1997. We 
have used LOTTOS considering 4x4 voxels of 6' in latitude and 7' in longitude, 
41 layers of 350 m height and 30-minutes batches for the Kalman filter. The 
tomographic fields have then been used to generate the reconstructed PWD 
which, in turn, have been processed to calculate the horizontal gradients and 
ZWD for each station. These have been compared to those obtained with GOA 
II. In Figure || we show the magnitude and direction of the 24-h mean value 
of the gradient as calculated using GOA II (yellow) and using the tomographic 
solution (green). In Figure |I| we show the time series of these magnitudes for 
a particular station (PGF3) and in Table 1 we show the correlation values over 
time for each station. This shows that the reconstructed 4D field is consistent 
with the data. Vertical distribution is shown in Figure | for a given latitude 
and time. We have computed the values of N(f, t) profiles with data from the 
European Center for Medium-Range Weather Forecasts (ECMWF) analysis and 



compared with our tomographic profiles, averaging them to meet the ECMWF 
maps resolution of 0.5 degrees. Results are presented for three different hours 
(06h, 12h and 18h) in Figure ^ showing good agreement. 

4. Conclusions 

We have successfully obtained tomographic 4D images of the pseudo-wet re- 
fractivity for a local dense network. The concept of slant pseudo-wet delays 
has been introduced. Simulated fields have been used to validate and tune our 
LOTTOS software package. Finally, data from the Kilauea network have been 
processed and the reconstructed 4D fields of the refractivity have been com- 
pared with on-site estimated horizontal gradients and ECMWF vertical profiles, 
yielding good correlation in both cases. Although more work is needed in this 
area, our results provide the proof-of-concept of the tomographic inversion of 
tropospheric GPS data. 
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Station 


ZWD 


North 


East 


AHUP 


0.98 


0.80 


0.83 


GOPM 


0.96 


0.79 


0.64 


KAEP 


0.98 


0.78 


0.71 


KTPM 


0.97 


0.57 


0.72 


MANE 


0.98 


0.90 


0.70 


MLPM 


0.99 


0.55 


0.88 


NUPM 


0.97 


0.62 


0.53 


PANU 


0.99 


0.53 


0.43 


PGF1 


0.98 


0.83 


0.83 


PGF2 


0.97 


0.74 


0.78 


PGF3 


0.97 


0.86 


0.86 


PGF5 


0.97 


0.85 


0.84 


PGF6 


0.97 


0.88 


0.89 


PULU 


0.98 


0.56 


0.41 


SAND 


0.99 


0.61 


0.82 


UWEV 


0.99 


0.69 


0.82 



Table 1: 24-h correlation factors for the ZWD 
and the Total Horizontal Gradients in both 
components, North and East, for each station. 
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Figure 1: Simulation Results: Error (in contours, labeled at 0.2 cm/Km steps) in the 
reconstruction as a function of longitude (in degrees East), for a constant latitude. 
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Figure 2: Impact of white noise in the tomographic solution. 
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Figure 3: 24-h Reconstructed horizontal gradients (green) compared to GOA II solu- 
tions (yellow). Latitude in y-axis and Longitude in x-axis. 
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Figure 4: Time series of the ZWD and Gradients as computed using GOA II (solid) 
and the tomographic fields (dashed) for PGF3. 
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Figure 5: Image of the refractivity field for a latitude of 19.28° N, February 1st, in the 
Hawaii network, having Longitude in x-axis (degrees East) and height (Km) in y-axis. 
Color scale is in cm/Km. 




Figure 6: Comparison of tomographic solution (solid line) with ECMWF (diamonds) 
for hours 06h (top), 12h (middle) and 18h (bottom). 



